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In this pedagogical article, we present a simple direct matrix method for analytically 
computing the Jacobian of nonlinear algebraic equations that arise from the discretization 
of nonlinear integro-differential equations. The method is based on a formulation of the 
discretized equations in vector form using only matrix-vector products and component- 
wise operations. By applying simple matrix-based differentiation rules, the matrix form 
of the analytical Jacobian can be calculated with little more difficulty than that required 
when computing derivatives in single-variable calculus. After describing the direct matrix 
method, we present numerical experiments demonstrating the computational performance 
of the method, discuss its connection to the Newton-Kantorovich method, and apply it 
to illustrative ID and 2D example problems. MATLAB code is provided to demonstrate 
the low code complexity required by the method. 

Keywords: analyticalJacobian; numerical methods; matrix calculus; Newton's method; integro-differential 
equations 



1. Introduction 

Many numerical methods for solving nonlinear integro-differential equations require 
computation of the Jacobian for the system of algebraic equations that arises when the 
continuous problem is discretized. For example, any Newton's method calculation requires 
computation of the Jacobian (exactly or approximately) during each Newton iteration [ 
[281 [29] . Unfortunately, calculation of the Jacobian can be a time-consuming and error- 
prone procedure for both the computer and the scientific programmer. 

In this pedagogical article, we present a simple direct matrix method for calculating 
analytical Jacobians of discretized, nonlinear integro-differential equations. The direct 
matrix method produces the Jacobian for the discretized equations directly in matrix 
form without requiring calculation of individual matrix elements. The essential idea is 
to first write the discretized, integro-differential equation explicitly in terms of discrete 
operators {e.g., differentiation and quadrature matrices [[HI [221 Hi]) and then use simple 
matrix-based differentiation rules to calculate the Jacobian directly [ [TOl [121 IIS [IS] . The 
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key observation underlying this approach is that there is a tremendous amount of structure 
in the nonhnear algebraic equations that arise from the discretization of nonlinear integro- 
differential equations. By taking advantage of this structure, the calculation of analytical 
Jacobians is reduced to nearly the level of complexity required to compute derivatives of 
scalar, single variable functions. 

The operator-based approach of expressing and analyzing discretized differential equa- 
tions has been used implicitly by the scientific computing community for decades, espe- 
cially in the context of the Newton-Kantorovich and related methods [El [3, [30]. However, 
a direct matrix approach seems to have been first formally described in one place by 
Chen who presented a collection of rules formulated in terms of specially defined matrix 
products [ [ini [121 [13] • In addition to using the method to solve nonlinear partial differ- 
ential equations [[121 [Ej; Chen used his formulation of the method to develop several 
interesting theoretical results {e.g., stability analysis of numerical methods for nonlinear 
time-dependent problems) based on the observation that when a nonlinear differential 
equation only has polynomial nonlinearities, there is a very close relationship between the 
discretized nonlinear differential equation and its Jacobian [ [TOl [11] . 

Mathematically, the present formulation of the direct matrix method is equivalent to 
Chen's approach. However, rather than introducing special matrix products, we rely solely 
on standard linear algebra operations augmented by component-wise operations [e.g., the 
Hadamard or Schur product [ 25j). In addition, we have chosen to use MATLAB notation 
in our formulation because of its prevalence in modern scientific computing. Working 
in MATLAB notation has the added benefit of making it almost trivial to translate the 
analytical calculations into working numerical cod^. 

Another feature of our formulation, which is also present to some extent in [ [10] , is 
the emphasis on the analogy between calculation of Jacobians for discretized, nonlinear 
integro-differential equations and calculation of derivatives for scalar functions of a single 
variable. To help strengthen the analogy with single-variable calculus, we organize the 
operations required to compute a Jacobian as a short list of simple differentiation rules. 

This article is organized as follows. In the remainder of this section, we compare 
the direct matrix method with several common methods for computing Jacobians. In 
Section[2l we present the direct matrix method, including a discussion of its computational 
performance and its relation to the Newton-Kantorovich method. Finally, in Section [31 
we apply the direct matrix method to two examples (one ID and one 2D) from the 
field of electrochemical transport. To demonstrate the low code complexity required by 
the direct matrix method, MATLAB code for the example problems is provided in the 
appendices. Throughout our discussion, we will focus solely on collocation methods where 
the continuous and discrete forms of the integro-differential equation have essentially the 
same structure. However, it is important to recognize that the direct matrix method can 
also be used for Galerkin methods by applying it directly to the weak-form of the problem. 

1.1. Comparison with Common Methods for Computing Jacobians 

One common approach for obtaining the Jacobian of a discretized, nonlinear integro- 
differential equation is to compute it numerically using finite differences of the grid func- 



^With today's powerful desktop and laptop computers, MATLAB is quite capable of handling moderate- 
sized production work. 
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tion or expansion coefficient values [ EHl EH]- Unfortunately, numerical computation of 
the Jacobian can be time consuming for some problems. Depending on the numerical 
method, it might be possible to reduce the computational cost of a numerical Jacobian 
by taking advantage of the sparsity pattern in the Jacobian [ HH [201 [29j , but determining 
the sparsity pattern can be complicated for nontrivial problems. 

As an example, consider the Poisson-Nernst-Planck equations for electrochemical trans- 



port [Him EE]: 

^ = V-(Vc+ + c+V0) (1) 

^ = V-(Vc_-c_V0) (2) 

eV20 = -(c+-c_), (3) 



where c± are cation and anion concentrations, respectively, and (j) is the electric potential, 
and e is a dimensionless physical parameter related to the dielectric constant of the elec- 
trolyte. The ffist two of these are the Nernst-Planck equations for ion transport and are 
simply conservation laws for cations and anions [[331 [36] . The last equation is the Poisson 
equation [[27], which provides closure for the Nernst-Planck equations. Note that in (jSj), 
the local charge density has been written in terms of the individual ions, which are the 
only source of charge density in many electrochemical systems. ([ID - ([2|) form a nonlinear 
parabolic system of partial differential equations, which suggests that we use an implicit 
time-stepping scheme to numerically solve the equations. This choice, however, requires 
that at time t„+i, we solve a nonlinear system of equations for c^^'^^^ which depends on 
an auxiliary variable which is in turn related to c|^^^'' through the Poisson equation. 
Note that in order to numerically compute the Jacobian for the resulting nonlinear system 
of equations for c^J!~^^\ we must solve the Poisson equation for each perturbation to the 
current iterate of Cj?^^"*. Therefore, for a pseudospectral discretization of ([T|) - Q using 
N grid points, numerically computing the Jacobian requires 0{N^) operations, which is 
much higher than the 0{N'^) elements in the Jacobiaro. 

Using an analytical Jacobian is one way to avoid the computational cost associated with 
numerical Jacobians. In principle, it is straightforward to derive the analytical Jacobian 
for the system of algebraic equations that arises when a nonlinear integro-differential 
equation is discretized. Index notation (also known as tensor notation) is perhaps the 
most common technique used to calculate analytical Jacobians. The basic idea behind 
the index notation method is to write the discretized form of the differential equation 
using index notation and then use tensor calculus to compute individual matrix elements 
in the Jacobian. For example, for a finite-difference or pseudospectral discretization, the 
discretized equations can be written in the form: 

Fj(Mi,M2, • • • ,MAf) = (4) 

for i = 1,2, ... ,N where Ui and are the value of the solution and the discretized 
differential equation at the i-th grid point (or more generally, the i-th collocation point). 



^While low-order discretization of the equations do not show this same disparity in the computation time 
and the number of elements (requiring 0{N^) operations for 0{N^) elements), they typically require 
many more grid points to produce an accurate solution. 
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Boundary conditions are included in this formulation by using the discretized boundary 
conditions (rather than the integro-differential equation) at grid points on the boundary 
(or immediately adjacent to if no grid points reside on the boundary The ij-th element 
of the Jacobian, J, for (jlj) is simply the partial derivative of Fi with respect to Uj: Jij = 
dFi/duj. While simple and straightforward, index notation suffers from the disadvantage 
of being somewhat tedious and error-prone. The main challenge in using index notation 
is keeping track of all of the indices when writing out and computing partial derivatives 
of the discrete equations. 

Automatic differentiation [ [2H [26] offers an important alternative when exact Jaco- 
bians are desired. Because it generates code for computing the Jacobian directly from the 
code used to evaluate the residual, automatic differentiation completely eliminates the 
possibility of human error when deriving the exact Jacobian and implementing it in code. 
Recent developments have made automatic differentiation available in several common 
programming languages (including MATLAB [[37]). While useful, automatic differentia- 
tion still takes some effort to use and may not always generate the most compact, efficient 
code. Active development in this area will certainly continue to improve the usability of 
automatic differentiation software and the performance of generated code. 

The direct matrix method has several advantages over the methods discussed in this 
section. First, the direct matrix method yields a more accurate Jacobian than finite 
differences and generally in less time (see Section [2^ . Second, because the method is 
based on simple differentiation rules, the calculation is straightforward and less prone to 
error than the index notation approach. The differentiation rules also make it easier to 
calculate the Jacobian for differential equations which depend on auxiliary variables, such 
as ([1]) - ([3]). From a programming perspective, calculation of the Jacobian directly in 
matrix form facilitates implementation of numerical methods for nonlinear problems in 
languages that have built-in support for matrix and vector operations {e.g., MATLAB 
and Fortran 95). Finally, having the Jacobian available in matrix form can be useful for 
analyzing properties of numerical methods [ [TU] . While it may be possible to to convert 
the element-wise representation of the Jacobian derived using index notation for simple 
problems, this step can be challenging for more complex problem^. 

2. The Direct Matrix Method 

There are two basic ideas underlying the direct matrix method for calculating analyti- 
cal Jacobians of discretized, integro-differential equations. First, rather than writing the 
discretized, integro-differential equations at each of the collocation points in terms of indi- 
vidual elements of the solution vector, we write the entire system of equations as a single 
vector equation expressed explicitly in terms of matrix-vector products and component- 
wise multiplication {e.g., Hadamard products). Second, the analytical Jacobian for the 
discretized system of equations is computed directly in matrix form by using simple differ- 

■^Care must be exercised when imposing boundary conditions, especially when using pseudospectral meth- 
ods [[23]. 

^Interestingly, the conversion from element-wise to matrix representation of the Jacobian often reveals 
the close relationship between the Jacobian for the discrete equations and the underlying structure of the 
original integro-differential equation. 
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entiation rules that are reminiscent of those used to compute derivatives in single- variable 
calculus. In this section, we develop both of these ideas in detail. Towards the end of the 
section, we comment on the computational performance of the direct matrix method and 
its relationship to the Newton-Kantorovich method [[8]. 

2.1. Matrix-vector Representation of Discretized Equations 

Writing the discretized, nonlinear integro-differential equation explicitly in terms of 
basic linear algebra and component-wise algebraic operations is the initial step of the 
direct matrix method. Because of the similarities in the structure between the discrete 
and continuous forms of the equations, the procedure is very straightforward. First, 
convert all differentiation and integration operators into their discrete analogues. Since 
both of these operations are linear, they become multiplication of vectors representing 
grid functions by differentiation and quadrature matrices, respectively: 

du ^ , , 

— ^ D*u (5) 
dx 

udx — > Q * u, (6) 

where the hat accent indicates a discretized field variable and D and Q are the differen- 
tiation and quadrature matrices associated with the choice of computational grid. 

Next, convert all point- wise algebraic operations and function evaluations in the con- 
tinuous equations to component- wise algebraic operations and function evaluations in the 
discrete equations. Some examples of the conversion process include: 

* v) (7) 

(8) 

(9) 
(10) 

In these examples, we have adopted the MATLAB convention of using .op to represent 
component-wise application of the op operation. Also, note that we have abused notation 
for component-wise function evaluations - f{u) represents the vector 

{f{u,)J{u2),...J{uN)) (11) 

not an arbitrary vector function of the entire solution vector u. Throughout our discus- 
sion, we will indicate component-wise and general functions of u by using lowercase and 
uppercase variables, respectively. 

2.1.1. Differential and Integral Operators in Multiple Space Dimensions 

It is important to emphasize that the matrix-vector representation is not restricted to 
scalar field equations or functions of a single variable. Handling vector equations is simple 
- vector equations may be treated as systems of equations and vector operations may be 
expressed in component- wise form. 

The construction of differential and integral operators for functions of multiple variables 
is slightly more complicated but straightforward. First, we represent grid functions as a 



dv 


-> u .* (D 


dx 


sm{u) - 


->■ sin(u) 


- 


U.A2 


exp(-u) - 


exp(-u). 
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(3,2) (3,2) (3,3) (3,4) 



(2,1) (2,2) (2,3) (2,4) 



(1,1) 



(1,2) 



(1,3) 



(1,4) 



X 



Figure 1. Example of a 2d grid with four grid points in the x-direction and three grid 
points in the y-direction. 



ID vector by selecting an ordering of the grid points. Then, we derive the differentia- 
tion and integration matrices associated with this choice of ordering. While there is no 
unique mapping from a multi-dimensional grid to a ID vector, it is important to choose 
the ordering of the grid points carefully because it directly affects the ease with which 
differentiation and integration matrices can be derived. 

For problems on logically rectangular computational domains, the computational grid 
may be constructed as a Cartesian product of one-dimensional grids [HI]. The result is a 
structured grid (possibly non-uniform depending on the discretization in each coordinate 
direction). To represent a grid function as a ID vector, the most natural way to flatten 
the grid is by using a lexicographic order for the grid indices. For example, on the small 
3 by 4 Cartesian grid in Figured], we could order the grid function, u, one row at a time: 

U = {Un, Mi2, Ul3, Ul4, U21, ^22, ^23, ^24, ^31, ^32, ^33, ^34)^ (12) 

where Uij is the value of u at {xi,yj). With this choice of ordering, discrete partial 
derivative operators are given as Kronecker product of the differentiation matrices with 
identity matrices [SIl- For the example in Figured], the partial differentiation matrices 
are 

= h®D4 (13) 

Dy = D^^h (14) 

where Dn and /„ are the Id differentiation and identity matrices of size n and ® denotes 
the Kronecker product. Note that care must be taken to ensure that the order of the 
Kronecker products is consistent with the ordering of the grid function vector. 

To illustrate these ideas, let us consider the discretized forms of the curl and divergence 
of a vector field, U . Once we have chosen an ordering of the grid points and derived the 
corresponding partial differentiation matrices D^, Dy, and D^, the discrete curl operator 
may be easily written as 

Dy*U:, - D^* Uy 

D,*a,- * f/, I , (15) 

*Uy- Dy*U^ 
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where Ux, Uy, and Uz are the components of U. The discrete divergence operator is also 
straightforward to derive: 



Dx*Ux+ Dy*Uy + Dz* Uz. 



(16) 



Boundary Conditions 

For problems in multiple space dimensions, it is often convenient to include discretized 
boundary conditions in a matrix-vector representation by first breaking the full differen- 
tiation and integration matrices into multiple components. Each component is defined by 
the grid points that it contributes to and the grid points it receives contributions from. 
For instance, it may be convenient to decompose a differentiation matrix into four mutu- 
ally exclusive components that: (1) use interior points and contribute to interior points, 
(2) use interior points and contribute to boundary points, (3) use boundary points and 
contribute to interior points, and (4) use boundary points and contribute to boundary 
points. Splitting the differentiation and integration matrices into separate components 
can be helpful when computing residuals and Jacobians for the different types of grid 
points in the computational domain. 

Deriving these components is straightforward using zeroth-order restriction and prolon- 
gation matrices [[9]. A zeroth-order restriction matrix is a matrix of zeros and ones which 
extracts a desired subset of elements from a vector. A zeroth-order prolongation matrix 
is also a matrix of zeros and ones, but it injects the elements of a restricted vector into a 
desired subset of the elements of a full-length vector. If {ii, 12, ■ ■ ■ , im} are the (flattened) 
indices of a subset of points from grid points, then the associated restriction matrix, R, 
would be an m X A^ matrix, R, with ones at the positions (1, ii), (2, ^2), • • • , (^, im)- The 
associated prolongation matrix, P, that injects a vector of length m into the positions 
{ii, ^2, • • • , im} of a vector of length A^ is simply the transpose of R: P = R^ . For example, 
the restriction and projection matrices for the two interior points of the grid in Figure [1] 
are given by 



To derive the differentiation matrix that uses values from grid points A and contributes 
to grid points B, we use the prolongation matrix for set A to expand the restricted grid 
function associated with A into a full-length vector and the restriction matrix for set B 
to compress derivative grid function to the set B: 



Restricted integration matrices are derived in exactly the same manner. 

2.2. Simple Differentiation Rules for Computing Exact Jacobians 

Once the continuous equations have been put in a discretized form that is expressed ex- 
plicitly in terms of matrix-vector products and component-wise operations, the analytical 
Jacobian for the discretized equations can be calculated by applying a few simple matrix- 
based differentiation rules. Because the differentiation rules are expressed completely in 
matrix form without any reference to individual elements in the Jacobian matrix, they 



R 



000001000000 
000000100000 



P = R^ 



(17) 



Da^b = Rb* D* Pa = Rb* D* R^. 
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allow us to compute the Jacobian directly in matrix form. In this section, we list these 
differentiation rules, which are essentially results from multivariate and matrix calculus 
applied specifically to the structure of discretized integro-differential equations. 

2.2.1. Matrix- Vector Product Rule 

The Jacobian of a matrix- vector product (which corresponds to a linear operator acting 
on a function in the continuous equations) is just the matrix itself: 

^iA*u)=A. (19) 

For example, the Jacobian of the discretized derivative of u, D * u, is just D. 

2.2.2. Diagonal Rule 

The Jacobian of a component-wise function / of a grid function u is a diagonal matrix 
with diagonal entries given by /'(«): 

^ = diag(m). (20) 

In essence, the diagonal rule is a way to use matrix notation to represent the fact that 
the differential in the i-th component of f{u) only depends on the change in the z-th 
component of u and is given by Sf{ui) = f'{ui)6ui, As an example, the Jacobian of sin({i) 
is diag (cos(u)). 

2.2.3. Chain Rules 

The Jacobian of a matrix A times an arbitrary function, F, of all of the components 
of M is A times the Jacobian of F: 

-(A,Fm-A,- (21) 

Similarly, the Jacobian of a function, F{u), when its argument is a matrix A times the 
grid function u is the Jacobian of F evaluated ai A*'u times A: 



§zFiA*u) 



9F , , 



* A (22) 



These rules are simply the chain rules for vector fields from multivariate calculus [[T]. 

For the special but common case when F is a component- wise function F{u) = f{u), 
(I2TI) reduces to A times the diagonal matrix with f'{u) on the diagonal: 

^{A*f{u)) = A*dmg{f'{u)) (23) 

and fl22l) reduces to the diagonal matrix with f'{A * u) on the diagonal times the matrix 
A: 

^/(A*M) = diag(f(A*M))*A (24) 
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2.2.4. Product Rule 

To compute the Jacobian of a component-wise product of general functions F and G 
of a grid function it, we use the product rule: 



— .* G(n))=diag(G(n))* — + diag(F(n))*— . 



(25) 



The derivation of the product rule follows from the expression for the variation of the i-th 
component of F{u)G{u) 



5{F,{u)G;{u)) 



Gi(u)^:rT-Su + Fi(u)^:rT-6u 



du du 



Gi{u)^:ri- + FAu, 



du 



du 



Su, (26) 



which yields the Jacobian 



Gi{u) (dFjdu) 
G2iu) {dF^/du) 

Gn{u) (dFi^/du) 



+ 



Fi{u) {dGi/du) 
F2{u) {dG2/du) 

Fn{u) (dGN/du) 



dW dG 
diag(G(^))* — +diag(F(«))*— (27) 



2.3. Example Jacobian Calculations 

As our first example, let us consider the ID Poisson equation: 



0. 



(2J 



To put this in discretized form, we need only replace the continuous second derivative 
operator by a discrete analogue: 



D'^*u 



P 



0. 



(29) 



Here, we have chosen to apply the discrete single derivative operator twice. Via a direct 
application of the matrix- vector product rule f[T^ . the Jacobian of the left-hand side 
of this equation is easily found to be D^. Since this is a linear equation, there would 
not normally be a need to compute the Jacobian of the left hand side of this equation. 
Moreover, the Jacobian for this example is very easy to calculate using alternative means 
(or even by inspection). We merely present it to illustrate the direct matrix method on a 
simple model problem. 

As a less trivial, let us calculate the Jacobian for the discretized form of the nonlinear 
function 



^2u 



du 
dx 



Converting this function to discrete form, we obtain 

f(u) = (e.A (2m)) .* {D*u). 
Using the product rule (!25l) . we find that the the Jacobian is given by 

d d(D * u) 

J = diag(D *u)* — (e.A (2u)) + diag(e.A (2n)) * --— — ^. 

du du 



(30) 



(31) 



(32) 
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Then applying the diagonal rule ( l20i) and the matrix- vector product rule ( |T9l) . we find 
that 

J = 2diag {D*u)* diag (e. A (2m)) + diag (e. A (2m)) * D, (33) 

which can be simplified to 

J = 2diag((D*M) .* (e.A (2m))) + diag(e.A (2m)) (34) 

by observing that diag (m .* v) = diag (m) * diag (m) . 

As a final example, let us calculate the Jacobian for the nonlinear algebraic equations 
that arise when solving the one-dimensional version of ([1]) - ([3]) using a simple backwards 
Euler discretization in time. Using the direct matrix approach for the spatial discretiza- 
tion, the nonlinear algebraic equations for c!^^^"* and cl'^^^ that need to be solved at each 
time step are: 

cf - At * a(^«+i) + D * (cf+') . * (D * 0)) ) - = (35) 

_ (^Jj2 ^ g{n+l) _ ^ ^ j^g(n+l) _ ^ ^ 0) j j _ ^L"^ = (36) 

eD%0+ (cV"+')-cL"+'^) = (37) 

where c|^^ are the concentrations at the current time step and At is the time step size. 
It is important to mention that several of the rows in ( 1371) will typically be replaced to 
impose the discretized form of the boundary conditions for 0. For illustrative purposes, 
let us suppose that we have simple Dirichlet boundary conditions for 0. In this situation, 
( 1371) is only imposed at interior grid points [ HI] . 

Using the simple differentiation rules from the previous section, the Jacobian of (!35ll 
with respect to c^^^"* is 

I-At(^D^ + D* diag (d*4>^+D* diag (|cf *D* j , (38) 

where / is the identity matrix and (^dc/) / dc^^^''^ is the Jacobian of with respect to c^^~^^\ 
To eliminate (d(f)/dc^_^^^^j from this expression, we simply apply the differentiation rules 



to (1571) with two rows eliminated for the boundary conditions and solve for the interior 
portion of (d(j)/dci^~^^^^: 

-(«^)":. (39) 



^^(n+l) / ^\ Jint 

+ / int 

where {D'^)int submatrix of that remains when all of the columns and rows 

corresponding to boundary grid points have been removed. Since the boundary values of 
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are fixed and the values of c^''"^'* at tlie boundaries do not affect the potential in the 
interior, the full Jacobian (^dc/) / dc^^^''^ is given by 



del) 



dc 



(n+l) 



■■■ 
■■• 



(40) 



where we have assumed that the first and last grid points correspond to boundary points. 

It is important to recognize that the form for ^c}0/9c^^^'' j in fl40p is specific to problems 

with Dirichlet boundary conditions for 0. For other boundary conditions, the inversion 

of the equation for / dc^^~^^^^ generally leads to different forms for the Jacobian. 

The Jacobian for (|35|) can now be explicitly computed by substituting fj40|) into fl38l) . 
The similar expression for the Jacobian of fl36l) is obtained using an analogous procedure. 
Using the direct matrix approach, we have reduced the calculation of the Jacobian to 
0{N^) (cost of matrix- inversion and matrix- matrix multiplies) compared to the 0(A^^) 
cost for computing a numerical Jacobian for high-order spatial discretizations. It is worth 
pointing out that in this example, the Jacobian for the concentrations does not depend 
explicitly on because the Poisson equation is linear. As a result, there is no need to 
solve for in order to compute the Jacobians for (l35l) and (1361) . For general problems, 
the Jacobian may depend on the auxiliary variable, so it might be necessary to solve the 
constraint equation. However, because only one solve for the auxiliary variables is required 
with the direct matrix method, the cost of computing the Jacobian is still dramatically 
reduced compared to using finite differences. 

2.4. Computational Performance 

In general, using the direct matrix method to compute a Jacobian is faster than cal- 
culating a numerical Jacobian. As mentioned in the previous section, the performance 
difference is expected to be large when auxiliary variables are involved in the expression 
of the residual. However, the direct matrix method yields higher performance even for 
problems where the residual is relatively simple. 

Figure [2] compares the performance of the direct matrix method against the MATLAB 
numjacO function for the two example problems discussed in Section [31 As we can see, 
the direct matrix method is at least an order of magnitude faster for both the ID and 
2D problems. For the 2D problem, the direct matrix method also shows superior scaling 
with the grid size. To ensure a fair comparison, we vectorized the residual calculation 
to minimize the number of function calls required by numjacO and avoided the use 
of matrix multiplication^, whenever possible, which benefited both methods. Matrix- 
matrix multiplications are especially detrimental for the direct matrix method because 
they can worsen the scaling of the Jacobian construction time with grid size to the point 
where the numerical Jacobian is faster to compute. For instance, in the left graph in 
Figure [2], a Jacobian computed using the direct matrix method with explicit matrix- 
matrix multiplications take 0{N^) time, which negates the performance benefits of the 



^For example, we express matrix- vector products of the form diag {u)*v as component-wise multiplication 
of two grid functions u .* v. 
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Figure 2. Comparison of the computational performance of direct matrix method (circles) 
and numerical Jacobian computed using MATLAB numjacO function (squares) for the ID 
electrochemical thin-film problem (Section 13. ip and the 2D metal colloid sphere problem 
(Section 13.21) . For both comparisons, the numjacO and direct matrix method codes 
were optimized by vectorizing the residual functions and avoiding matrix multiplications 
whenever possible. The data for these graphs were generated on a 2.4 GHz MacBook Pro. 



method compared with a numjacO implementation before even reaches 1000. 

In addition to avoiding matrix-matrix multiplication, it is important to use sparse ma- 
trices when possible. For problems in more than one space dimension, sparse matrices are 
produced when Kronecker products with identity matrices are used to construct differenti- 
ation matrices even if the ID differentiation matrices are dense. Not only does the memory 
required for dense matrix representations easily exhaust the memory on workstations and 
laptops, dense matrix representations also leads to poor computational performance when 
applying and multiplying the matrices. In general, sparse matrix operations have better 
scaling properties as the grid size grows. 

2.5. Relationship to the Newton-Kantorovich Method 

The direct matrix method for computing the Jacobian of discretized integro-differential 
equations is closely related to the calculation of the Frechet derivativ^ used in the Newton- 
Kantorovich method [[8l[30] (also known as quasilinearization [[21]). The basic idea behind 
solving nonlinear integro-differential equations using the Newton-Kantorovich method is 
to carry out Newton's method in function space. For each Newton iteration, we compute 
the Frechet derivative of the integro-differential equation in function space and numerically 
solve the resulting linear integro-differential equation for the correction to the current 
iterate of the solution. Essentially, the Newton-Kantorovich method reverses the order 
of (1) discretization of the continuous problem and (2) Newton iteration. Because the 
equations to be solved during each Newton iteration is linear, there is no need to compute 

^Recall that the Frechet derivative for nonhnear functionals is the generahzation of the Jacobian for 
nonUnear functions over finite-dimensional spaces [ [8j [STJ |40] . For intuition, Ortega and Rheinboldt 
provide a nice discussion of Frechet derivatives in the context of finite-dimensional spaces [ [M] . 
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a Jacobian of tlie discretized equations. 

An important feature of the Newton-Kantorovich method is that the numerical dis- 
cretization used to solve the linearized equation during each Newton iteration can, in 
principle, be completely independent of the discretization used to compute the residual of 
the nonlinear integro-differential equation. This freedom can affect the convergence rate 
of the method depending on the degree to which the discretized form of the linearized 
problem approximates the Jacobian of the discretized residual equation. 

Because the direct matrix method begins with a discrete equation possessing the same 
mathematical structure as the continuous residual equation, it produces a Jacobian that 
is a discrete analogue of the Frechet derivative for the continuous integro-differential 
equation. Unlike the Newton-Kantorovich method, however, the direct matrix method 
produces the unique Jacobian associated with the particular choice of discretization for 
the residual of the nonlinear integro-differential equation. The freedom to independently 
choose the numerical discretizations for the residual equation and the Frechet derivative is 
not present in the direct matrix method. As a result, given a numerical discretization for 
the residual equation, the direct matrix method can be viewed as a way to generate the 
optimal discretization for the linearized equation that arises during each Newton iteration 
of the Newton-Kantorovich method. 



3. Applications 

Analysis of electrochemical systems is a classical subject that has recently seen a renewal 
of interest. Modern electrochemical systems of interest include ion channels in biological 
membranes [[2l[3l[35], microfluidic devices based on electro-osmotic flows [El [38], and thin- 
film battery technologies [ [32l [391 US] • A common feature of many of these applications 
is that the electrochemical system is operated under extreme conditions, such as large 
applied fields or very small physical size [[H [IS]- In these regimes, numerical solutions of 
the nonlinear governing equations are useful for gaining insight into the rich behavior of 
these systems. As we shall see, the direct matrix method makes it easy to compute the 
analytical Jacobian required to solve these nonlinear equations using Newton's method. 

3.1. Electrochemical Thin-Films 

Analysis of ID electrochemical systems leads to an example of a nonlinear integro- 
differential equation. For steady-state electrochemical thin-films made up of a dilute 
solution of symmetric binary electrolyte with faradaic reactions at the surfaces of the 
thin-film [HI US], the electric field, E, satisfies the second-order differential equatior0 



^This equation is mathematically equivalent to the Poisson-Nernst-Planck equation formulation of elec- 
trochemical transport [[4]. To simplify the discussion, equation ((4T|) is a slightly modified form of the 
master equation in [HI [12] derived by making the substitutions x ^ {x + l)/2 and E 2E. 



14 



K. T. Chu 



on the domain (—1, 1) subject to boundary conditions that represent the kinetics of elec- 
trode reactions 

-fc,(c(l)+p(l))+j,-j = (42) 
A;e(c(-l) + p(-l))-j,-j = 0, (43) 

where c is the average ion concentration, p is the charge density, j is the current density 
flowing through the thin- film, e is a parameter related to the dielectric constant, kc and 
jr are reaction rate constants, and cq is the following expression 



co = (l-j)+e' 



2E(1) -2E(-1) - / E^dx 



(44) 



The average ion concentration and charge density are related to the electric field via the 
equations 

dE 

c{x) = co+3{x + l) + 2e^E^ , p{x) = 4e'^. (45) 

We can solve this set of equations via Newton's method using a systematic application 
of the direct matrix method. To discretize the equations, we use a pseudospectral method 
based on the Chebyshev grid on the interval [—1, 1]. The differentiation matrix, D, for 
this computational grid is just the standard differentiation matrix for the Chebyshev 
grid [ El [221 SI]- Foi' numerical integration, we use the Clenshaw-Curtis quadrature 
weights [E], which we denote by the row vector w. The quadrature weights are used 
to construct a quadrature matrix, Q, which is the analogue of the differentiation matrix: 
Q = [w'^, , . . . , w'^Y . When a grid function / is multiplied by Q, the result is a vector 
where all entries are equal to the numerical approximation of the integral of /. 

With these discrete operators, we can put (HTj) into matrix-vector form: 



with 

Co = (1 - j) + (2Ei -2Em-Q* {E.A 2)) , (47) 

where we have chosen to order the indices so that xi = 1 and xa? = — 1 (this follows the 
convention used in [UT] and in the code in Appendix |A]) . The boundary conditions are 
imposed by replacing the discrete equations corresponding to xi and xn with 

- h (co + 2j + (2El + 4Di * ) + j,. -j = (48) 
K f Co + (2EJ, + AD,,*e))- jV - j = 0, (49) 



{D'*E--E. A3 \ - - [Co + 3{x + 1)) . * E - = (46) 



where Di and are the rows of the differentiation matrix corresponding to Xi and xn, 
respectively, and Cq is a single component of Cq. 
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Figure 3. Solution of electrochemical thin-film equations (HTj) - (H3|) computed using 100 
grid points with e = 0.01, kc = 10, and jr = 10 for j = 1.5 (solid), j = 1.0 (dash), and 
j = 0.5 (dot-dash). 



The Jacobians for the left-hand side of these discrete equations are now easily computed. 
Applying the differentiation rules from Section [2l the Jacobian for the interior grid points 
is 



Jint = e' (D' - -dmg[E.A2 



-diag (Co + j(x+ i; 



-diae i E] * 



dCo 
dE 



(50) 



with 



dE 



VL 



2 ■■■ 



2 



-2 



- 2Q * diag E 



J 



(51) 



The Jacobian for the discretized boundary conditions are similarly calculated: 

'dco 



Ji 

Jn 



-k 



dE 



+ Ae'[Ei ... 0] +4e^L>i 



^ + 4e2[0 ... QEN]+Ae^D 
dE 



N 



where 



dE 



eM [2 ... -2]-2w* diag E 



(52) 
(53) 

(54) 



From the perspective of computational performance, the above formulation of the Ja- 
cobian is suboptimal because it includes a matrix-matrix multiply in fl5U]) that can be 
avoided. To reduce the time required to compute the Jacobian, the key observation is 

that each row of ^ is equal ^. Therefore, diag [e]*^ is more efficiently computed as 



dE 



dE 
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Figure 4. Plots of the absolute value of the spectral coefficients, a„, for the numerical 
solution of the electrochemical thin-film equations with j = 1.5 (solid), j = 1.0 (dash), 
and j = 0.5 (dot-dash) as a function of the basis function degree when ( HTi) - ( 143|) are 
solved directly (left) and using the variable transformation (155|) to place more grid points 
in the boundary layers (right). Although the convergence rate as a function of number 
of grid points is geometric in all cases, the variable transformation significantly reduces 
the number of grid points required obtain the best solution possible given the precision 
of the computation. Note that the spectral coefficients were calculated from a numerical 
solution generated using 250 grid points to make the roundoff plateau more apparent. 



the Kronecker product of E and The evaluation of the residual can also be improved 

by recognizing that all of the elements of Co are equal to cq, but this optimization has a 
far smaller impact than the reformulation of the Jacobian. 

Now that we have explicitly computed all of the components for Newton's method, it 
is straightforward to write a program to solve the nonlinear integro-differential equations 
for electrochemical thin-films. The MATLAB code for solving is relatively short and 
runs quickly (see Appendix |A]). One special issue that arises for this problem is that 
continuation methods [ [8] are required to obtain good initial iterates for the Newton 
iteration at high current densities. Figure [3] shows the numerical solution of fl4ip - (1431) 
computed using 100 grid points with e = 0.01, kc = 10, and jr = 10 for various values 
of j. As expected, we observe geometric convergence with respect to the number of grid 
points (see Figure H]). Notice that at higher current densities, we see slower convergence 
rates due to the presence of greater structure in the solution. 

While quite satisfactory, the convergence rate for the numerical discretization - 
(jinD as a function of the number of grid points is limited by the need to resolve the 
boundary layers. By using a mapping of the computational domain that allows us to 
place a few grid points within the boundary layers, we can obtain a faster convergence 
rate. For example, by using the variable transformation: 

X = — tanh (atanh(/?)?/) , E{x) = ^ cosh^ (atanh(/?)|/) E{y), (55) 

where f3 is an adjustable parameter less than 1, we can significantly reduce the number 
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of grid points required obtain a solution accurate to machine precision (see Figure HI). 
It is interesting to note that the optimal value for (3 depends on the current density j. 
For j = 0.5 and j = 1.0, a P value of 0.9 yields near optimal results. For j = 1.5, 
however, the fastest convergence is obtained near (3 = 0.75. As is typical, the transformed 
electrochemical thin-film equations are bit more complicated to deal with than the original 
equations. However, the direct matrix method makes it straightforward to discretize 
the transformed equations and compute the exact Jacobian for the resulting nonlinear 
algebraic equations (see Appendix [B|. 

3.2. Double Layer Charging of Metal Colloid Sphere at High Applied Electric 
Fields 

Analysis of double layer charging for colloid systems subject to applied electric fields 
gives rise to nonlinear differential equations in multiple space dimensions with complicated 
boundary conditions. In the electroneutral limit [HI [331 ES], the steady-state governing 
equations for systems composed of symmetric binary electrolyte are [ [T7] 



V^c = (56) 

V-(cV0) = 0, (57) 

where c is the average ion concentration and (p is the electric potential. For metal colloid 
surfaces, the appropriate boundary conditions are [ [TTJ [18] 

= eVs ■ (gV,lnc + u;V,</)) -c^ (58) 

on 

Oc 

= eV, ■ (wV,lnc + gV,0) - — (59) 

on 

q = -2v/csinh(C/2) (60) 

w = 4v^sinh2(C/4) (61) 

v-(P = C + 25ycsinh(C/2) (62) 



where q and w are the excess charge and ion concentration in the boundary layer, ( is the 
electric potential drop across the boundary layer, v is the potential of the metal colloid, 
and ^ is a parameter related to the capacitance of the boundary layer. 

As a model problem, we solve these equations for a metal colloid sphere subjected to 
a uniform applied electric field of strength E in the ^-direction. To avoid infinite values 
of the electric potential, the numerical model is formulated in terms of = (p + Ez, the 
deviation of the electric potential from that of the uniform applied field. The spherical 
geometry of the problem also allows us to demonstrate the use of the direct matrix method 
on a non-Cartesian (though still logically rectangular) grid. 

While this problem may seem daunting, it is straightforward to obtain a solution numer- 
ically by using Newton's method with an analytical Jacobian computed using the direct 
matrix method. Taking advantage of azimuthal symmetry, we discretize the equations in 
spherical coordinates on a 2D pseudospectral grid that is the tensor product of grids in 
the radial and polar angle directions. We use a shifted semi-infinite rational Chebyshev 
grid [ [8] in the radial direction and a uniformly spaced grid for the polar angle direction. 
The required differentiation matrices are constructed using Kronecker products, and the 
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boundary conditions are handled using restriction and prolongation matrices as discussed 
in Section [2.1.1[ 

To facilitate the formulation of the matrix- vector representation of the equations, let us 
fix our notation. Let Dj. and Dg be the radial and angular contributions to the discrete 
divergence operator, Gr and Gq be the radial and angular components of the discrete 
gradient operator, and L be the discrete Laplacian operator. Also, let n and s subscripts 
denote normal and tangential derivative operators at the surface of the sphere. 

For the purpose of discussion (and implementation), it is convenient to decompose the 
discrete differential operators into pieces that correspond to contributions from finite and 
infinite grid points. For example, L can be decomposed into L-^ and L°° which respectively 
account for the contributions to the Laplacian operator from finite and infinite grid points; 
that is, L * c = * Cf + L°° * Coo, where c/ and Coo are the concentration values at finite 
and infinite grid points respectively. Similarly, to impose the boundary conditions, we use 
derivative operators that act only on surface values. Surface operators and surface field 
values will be denoted with superscripts s and subscripts s, respectively. Finally, to refer 
to values at interior grid points {i.e., finite grid points that are not on the surface of the 
sphere), we use the subscript i. 

In this notation, the discretized form of the bulk equations dSHD and (|Fr|) are given by 



= Fi = * c/ + L~ * Coo 
= F2 = Dl* 

+ dI* 



Cf .* yG{. * ipf — E cos 6 
Cf .* (Gl*'npf + Esm9 



-D^ * (coo .* EcosO) 
* (coo .* EsinO) . 



(63) 
(64) 



In these equations, the unknowns are the values of the c and at finite grid points; values 
at infinity are specified by the boundary conditions and so are known quantities (which 
is why Coo does not have a hat accent and tpoc, = does not show up at all). In discretized 
form, the boundary conditions on the surface of the sphere are 



= Hi 



= Ho 



eDs* 



q .* {G"" *\ncs) + w .* (^G''*ips 
Cs .* (^GI* ipf + E cos9^ 

w .* (G"" * Incs) + g .* (G''*iIjs 



eDs* 



G'*Ecos6 



G'*Ecos6 



(G^ c^ + * Coo) . 



(65) 
(66) 



Closure for these equations is given by using fIBU]) - fl^ to relate q and w to the zeta- 
potential and using (162|) to compute the zeta-potential from and Cg. 

The direct matrix method makes it straightforward to derive the analytical Jacobian 
for the system of equations (I6OII - (l66ll . The derivatives of Fi and F2 with respect to the 
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DFi 
DFi 

DF2 

Dcf 
DF2 



dI * diag (gI * ipf - EcosO^ + * diag (gI *ipf + EsinO 
Dl * diag (c/) * Gl + * diag (cj) * 



(67) 
(68) 

(69) 
(70) 



The derivatives for the discretized boundary conditions are more complicated because g, 
w, and Cs implicitly depend on the unknown variables and because surface grid points 
must be treated differently than interior grid points. However, a systematic application 
of the differentiation rules in Section 12.21 yields the analytical Jacobian directly in matrix 
form: 



DHi 



DHi 
DHi 



= -Ds * diag {q ./ Cs .* (G"" * In c^) ) 

/ dC 

- eD, * diag ( .* cosh(C/2) .* — .* (G'*lnc,) 

+ eDs * diag(g) * G"^ * diag (1 ./ c^) 

+ -Ds * diag{ w ./ Cs .* {G"" * ips - G" * E cos 6) ) 

+ eDs * diag ( JFs .* sinh(C/2) .* .* (G' * ips - G' * E cos 6) 
V 9cs 

- diag {Gl *4jf + E cos 9) 
= 



(71) 
(72) 



—eZ^s * diag -v/cl •* cosh(C/2) .* tt— •* (G"**lnCs) + eZ^^ * diag(w) * G'^ 

dtps 



DHi 



+ eD, *diag ( .* sinh (C/2) .* ^ .* iG' * tps - G' * E cosO) 
- diag {cs) * Gl 
= -diag(c,)Gt. 



(73) 
(74) 
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DH2 



= -D,* diag{w ./ c, .* (G^^lnc^) ) 

+ eD, *diag(^v^.* sinh(C/2) .* ^ .* *\nCs] 

+ eDs * diag(ty) * G"* * diag (1 ./ Cg) 

+ ^Ds*dmg{q ./ Cs .* {G' *iJs~G' *E cos 9) ) 

/ dC 
- eDs * diag { y/c^ .* cosh(C/2) .* t;— ■* {G'' * tps - G"" * E cos i 

\ OCg 



DH2 

Da 
DH2 



DH2 

D^, 

where 



-Gi 



= eDs * dia.g \^^/c^ .* sinh(C/2) .* 
— eDs * diag [ y/c^ .* cosh(C/2) .* 



dtps 
d^jjs 



.* {G'iJs-G'*Ecos9) 



dips 

K 

dcs 



1 



l + (5V^cosh (C/2) 
(5sinh(C/2) 



(75) 
(76) 

* {G' * Incs) ) +eDs* diag(g) * G' 



(77) 
(78) 

(79) 
(80) 



TcTfl +5v^cosh (C/2)]' 
The Jacobian for the system of equations is obtained by assembhng these pieces: 



J 



r dFi dF\ 

dc 

dF2 dh 
dc 

dHi dHi 
dc 

dH2 dH2 
L dc dip 



where the Jacobians for Hi and H2 are constructed from 
operators. For instance, 



^1) 



I]) - fl78|) using restriction 



dHi 
dc 



dHi 



dHi 

OCi 



(82) 



where Rs and Ri are restriction operators for surface and interior grid points, respectively. 
While the formulas may look complicated to program, they are actually quite easy to 
implement in MATLAB (see Appendix lUj) . 

Figure O shows numerical solutions obtained using the above residual and Jacobian 
formulas. As for the electrochemical thin-film example, continuation is required to obtain 
good initial iterates for the Newton iteration at high values of the applied electric fields. 
The solutions shown are computed for E = 10, v = 0, e = 0.01, and 6 = 1 using 30 
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Figure 5. Solution of equations flS^ - (jHSD computed using 30 grid points in both the 
radial and polar angle directions with E = 10, v = 0, e = 0.01, and S = 1. 



grid points in both the radial and polar angle directions with scale parameter set to 0.5 
for the shifted rational Chebyshev grid. Using pseudospectral grids and the analytical 
Jacobian, the solution is obtained very quickly, requiring only a few Newton iterations 
for each continuation stage (and less than a minute of computation time on a 2.4 GHz 
MacBook Pro). 

4. Conclusions 

In this article, we have presented a direct matrix method for calculating analytical 
Jacobians for discretized, nonlinear integro-differential equations. Because this method is 
based on simple matrix-based differentiation rules, it is less tedious and less error prone 
than other approaches for computing analytical Jacobians. Furthermore, because it yields 
the Jacobian in matrix form, it is very easy to use languages that support vectorized 
computation to implement numerical methods that require the Jacobian. 

One interesting possibility that the direct matrix method presents is development of 
high-level automatic differentiation tool for discretized nonlinear integro-differential equa- 
tions. In contrast to traditional automatic differentiation methods [I211I2S] which operate 
at the level of individual scalar operations, automatic differentiation methods based on 
the direct matrix method would operate on the discrete differential operators associated 
with the continuous differential equation. Such an automatic differentiation tool could 
be useful for completely eliminating the need for a researcher to compute the Jacobian of 
discretized nonlinear integro-differential equations by hand. 
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A. MATLAB Code for Electrochemical Thin-Film Example 



This code relies on cheb.m and clencurt.m [ST]. 

% parameters 

mysmmmmmmmmmmmmmmmmmmmm^ 

N ^ 200; 

j ^ 1.5; cpsilon ^ 0.01; k_c ^ 10; j_r ^ 10; 
res_tol — Ic— 8; inax_iters — 20; 



% compute grid , differentiation matrix , and quadrature weights 

[D,x] — cheb(N— 1); % Chebyshev differentiation matrix 

L — D*D; %Laplacian operator 

[x,w] — clcncurt(N— 1); % Clcnshaw — Curtis quadrature weights 



% set up continuation in j 

--/fi-^-/ri-fn-frr-fri--/ri-/ri-/ri-fri-fri--/ri-^-/ri-f^ ■ 

j_start — 0.5; dj — 0.1; j_cur — j_start; 



% generate initial iterate for Newton iteration 



cO = 1— j_cur; c — cO + j_cur*(x+l); E— — 2*j_cur./(j_cur*(x + l) + c0); 



% Newton iteration with simple continuation 

while ( j.cur <= j & dj > ) 

% display j _cur 
j_cur — j_cur 

% initialize Newton iteration 

count — 0; % iteration count for single Newton iteration 

cO = 1-j.cur + opsilon - 2*{2*E(1) -2*E(N)-w* (E . - 2 ) ) ; 

res = cpsilon " 2* (L*E — 0.5*E . " 3 ) ... 

- 0. 25 * ( cO + j -cur • (x + 1)) . *E - 0.25»j.cur; 

res{l) = -k_c*(cO + 2*j_cur + epsilon"2*{2«E(l)"2 + 4*D(l ,:)*E)) 

+ j_r — j_cur; 

res (N) = k_c * (cO+epsilon - 2* (2*E(N) -2 + 4*D(N, : ) * E) ) ... 

— j -r - j _cur ; 
res_norm — norm( res , inf ) ; 

while ( {res_norm > res_tol) & (count < max_iters) ) 

% construct Jacobian for interior grid points 
dcO.dE = — 2* c psi Ion "2»w. *E ' ; 
dcO.dE(l) = dcO.dE(l) + 2*cpsilon"2; 
dcO.dE(N) = dcO.dE(N) - 2* cpsilon "2; 
J — epsilon"2*L ... 

- diag { 1 . 5* epsilon - 2» (E. *E) + . 2 5 * ( cO+j .c u r * ( x + 1 ) ) ) ... 

- 0. 25* kron (E, dcO.dE ) ; 

% construct Jacobian for boundary conditions 
J(l ,:) = -k_c*(dc0.dE + 4* c p s i 1 o n ■ 2 *D ( 1 ,;)); 
J(l,l) = J(l,l) - 4* k_c* epsilon ■2*E( 1 ) ; 
J(N,:) = k_c*(dc0.dE + 4* e p s i 1 o n " 2 *D (N , : ) ) I 
J(N,N) = .I{N,N) + 4* k_c* epsilon ■2*E(N) ; 



% compute delta_E and update solution 
delta.E = -.I\res; E = E + delta.E; 



% update residual 

cO = 1-j.cur + epsilon - 2* {2*E(1) -2*E(N)-w*{E. - 2 ) ) ; 

res = epsilon " 2* (L*E-0.5*E. " 3) ... 

~ 0. 25 * ( cO+j _cur * (x + 1)) . *E - 0.25*j.cur; 
res{l) =-k_c * ( cO +2* j _ c u r + c p s i 1 o n ■2*(2*E(1)'2 + 4*D(1 ,:)*E)) 

+ j_r — j_cur; 
res{N) = k_c * ( cO+epsilon - 2* (2*E(N) -2 + 4*D(N, ; ) *E) ) ... 

- j -r — j _cur ; 

% update loop variables 
res_norm — norm(res , inf) 
count — count + 1 

end % Newton iteration loop 



% update continuation variables 
if (j - j-cur < dj ) 

dj = j - j-cur; 
end 

j_cur — j_cur + dj ; 



% plot solution 

figure(l); clf; 
plot{x,E,'k— '); 
axis([-l 1 -100 0]); 

xlabcl('x'): ylabel('E' , 'Rotation' ,0); 
% plot spectral coefficients 

coefs = abs{fft([E; f 1 i p u d { E { 2 ; end - 1 ) ) ] ) ) ; 
figure (2) ; clf ; 
scmilogy(cocfs (1:N) . 'ko'); 
axis ( [0 N le-15 le4 ] ) ; 

xlabel('n'); ylabel('|a_n| ', 'Rotation' ,0,' Position ',[—32 5e— 6]); 



B. MATLAB Code for Electrochemical Thin-Film Example 
with Variable Transformation 



This code relies on cheb.m and clencurt.m [141]. 



% parameters 

N — 200; beta — 0.75; alpha — atanh(bcta); 
j ^ 1.5; epsilon ^ 0.01; k_c ^ 10; j _r ^ 10; 
res_tol — le— 8; max_iters — 20; 



% compute grid , differentiation matrix , and quadrature weights 

[D,y] — cheb(N— 1); % Chebyshev differentiation matrix 

[Yiw] — clcncurt(N— 1); % Clenshaw — Curtis quadrature weights 

X— tanh(alpha'(;y)/bcta; % mapped grid points 

gamma — bcta/alpha*cosh{alpha=i=y).''2; % derivative transformation factor 

diag_gamma — d i ag ( gamma ); %cachediag( gamma) 

L — diag_gamma *D* diag_gamma *D ; % trans formed Laplacian operator 



% set up continuation in j 

%rsmmmmmmmmmmmm 

j_start — 0.5; dj — 0.1; j_cur — j_start; 



y^e/ir/(r/<v&ir/(r/<v&ir/(r/<v^ 

% generate initial iterate for Newton iteration 

mysmmmmmmmmmmmmmmmmmmsmms^^ 

cO — 1— j_cur; c — cO + j_cur*(x+l); E — — 2* j _ c u r . / ( j _c u r * ( x + l) + cO ) ; 



% Newton iteration with simple continuation 

WMmmmmmmmmmmmmm 

while ( j_cur <- j & dj > ) 

% display j_cur 
j_cur — j_cur 

% initialize Newton iteration 

count — 0; % iteration count for single Newton iteration 
cO — 1— j_cur ... 

4- epsilon " 2 * ( 2 * gamma ( 1 ) *E(1) — 2 *gamma{N) * E (N) -w* (gamma. * (E . ' 2 ) ) ) ; 
res — epsilon''2*(L*( gamma .^t^E) — 0.5*( gamma ."3).*(E."3)) ... 

— 0.25* ( cO + j _cur * (x + 1) ) gamma . *E— 0.25*j_cur; 
res(l) — — k_c*(c0 + 2*j_cur ... 

+ cpsilon ' 2* (2*gamma( 1)'2*E(1)'2 ... 

+4*gamma{ 1 ) *D ( 1 , : ) * ( gamma. * E ) ) ) ... 

+ j -r - j -c u r ; 

res (N) ^ k_c * ( cO + ep sil on " 2* (2*gamma{N) "2*E(N)"2 ... 

+4*gamma{N) *D(N , : ) * ( gamma. * E ) ) ) ... 

- j _r - j _cu r ; 
res_norm — norm{rcs,inf); 

while ( {rcs_norm > rcs_tol) &c (count < max_itcrs) ) 

% construct Jacobian for interior grid points 
d c _d E — — 2*epsilon"2*w.*( gamma . * E ) ' ; 
dcO_dE(l) ^ dcO_dE(l) + 2 * e p s i 1 o n " 2 * gamma ( 1 ) ; 
dcO_dE(N) ^ dcO_dE(N) - 2 * e p s i 1 o n ' 2 * gamma(N ) ; 
J — epsilon "2*L* diag_gamma . . . 



— diag(l. 5* epsilon ~2*( gamma . 3) . * (E. " 2) . . . rf^ 

+ 0.2 5*gamma. * ( cO+j _cur*(x + l))) ... 

— 0.25*kron (gamma. *E, dcO_dE ) ; 



% construct Jacobian for boundary conditions 

J(l ,:) — — k_c*(dcO_dE + 4* e p s i 1 o n ' 2 gamma ( 1 ) * D ( 1 , : ) * diag_gamma ) ; 
J(l,l} ^ J(l,l) - 4* k_c* epsilon "2*gamma( 1) 2*E ( 1 ) ; 

J(N,:) — k_c*(dcO_dE + 4* c p s i 1 o n " 2 * gamma (N) *D ( N , : ) * di ag.gamma ) ; 
J(N,N) ^ J(N,N) + 4*k_c*cpsilon "2*gamma(N)"2*E(N) ; 

% compute delta_E and update solution 
dclta.E ^ -,T\rcs; E ^ E + dclta.E; 

% update residual 
cO — 1— j _c u r ... 

+ epsilon ' 2 * ( 2 * gamma ( 1 ) *E(1) — 2*gamma(N) *E (N) -w* ( gamma . * (E . ' 2 ) ) ) ; 
res — epsilon''2*(L*( gamma .*E) — 0.5*( gamma .*3).*(E."3)) ... 

~ . 2 5 * ( cO+j _cur=i=(x + l)).* gamma. *E- 0.25*j_cur; 
res(l) — — k_c*(c0+2*j_cur ... 

+ epsilon " 2 * ( 2 * gamma ( 1 } "2*E(l}"2 ... 

+4*gamma( 1)*D(1 ,:)*( gamma . * E ) ) ) ... 

+ j-T — j-cur; 
res(N) ^ k_c * ( cO+epsilon " 2* (2*gamma(N) " 2*E(N) " 2 ... 

+4*gamma(N) *D(N, : ) * ( gamma. *E) ) ) ... 

— j -r — j -cur ; 

% update loop variables 
res_norm — norm(res , inf) 
count — count + 1 

end % Newton iteration loop 

% update continuation variables 
if (j - j-cur < dj ) 

dj ^ j - j_cur ; 
end 

j_cur ^ j-cur + dj ; 



end 



% plot solution 

figure (1); clf; 
plot(x,E.* gamma , ' k — ' ) ; 
axis ([-1 1 -100 0] ) ; 

xlabel ( 'x' ) ; ylabel ( 'E' , 'Rotation ' ,0); 
% plot spectral coefficients 

coefs = abs(fft([E; f 1 i p u d { E { 2 : end - 1 ) ) ] ) ) ; 

figure (2) ; clf ; 

semilogy( coefs (1:N) , 'ko '); 

axis ( [0 N lc-15 le4 ] ) ; 

xlabcl('n'); ylabel ( '|a_n| ', 'Rotation ' ,0, 'Position ',[—25 5e— 6]); 



C. MATLAB Code for Double Layer Charging of Metal Colloid 
Sphere 



This code relies on cheb.m [I41j. 



% Parameters 

% physical parameters 

V ^ 0; E ^ 10: epsilon ^ 0.01; delta ^ 1; 
% grid parameters 

N_r — 3 ; % number of grid points in radial direction 

N_t — 3 ; % number of grid points in polar angle direction 

L_r — 0.5; % scale parameter in radial direction 

% continuation parameters 
E_start ^ 1; E.final ^ E; dE ^ 0.5; 

% Newton iteration parameters 

res_tol — le— 8; delta_tol — le— 13; max_iters — 20; 
% zeta — potential iteration parameters 

zeta_max_iters — 20; zeta_delta_tol — le — 13; zeta_res_tol — le— 9 

% boundary conditions 
c_infinity — 1; 



% Construct computational grid and differentiation operators 



% construct the differentiation matrix for the radial coordinate 
[D_y , y ] = cheb (N_r ) ; 

one _mi nus_y = spdiags(l— y,0,N_r + l,N_r + l); 

D_r = 0.5/L_r*(one_min us_y " 2 ) * D_y ; 

warning off MATLAB: divideByZero 

r - L_r*(l + y)./(l-y); 

warning on AdATLAB : dividcByZcro 

r ^ r + 1; % shift to 1 

% construct the differentiation matrix for the polar angle coordinate 
theta ^ (2*[l:N_t]'-l)*pi/2/N_t; repmat(theta,l.N_t); 
c ^ ones { 1 , N_t ) . * ( - 1 ) . " { [1 : N_t ] + 1 ) ; 
off_diag_D ^ ... 

repmat(c,N_t ,1).* sin(N_t*T).*sin(T')./{ cos{T')-cos{T)+eye{N_t)); 
d iag_D — — 0.5*cot(theta); 

D.theta ^ t r i u { o f f _d i a g _D . 1 ) + t r i 1 ( o f f _d i ag _D , - 1 ) + d i ag ( d i ag_D ) ; 

num_gridpts_r — length(r); num _gr i dpts_theta — length(theta); 
num_gridpts — {num_gridpts_r— l)*N_t; 
num_gridpts_interior — num_gridpts— N_t; 

% cache common expressions 

one_over_r — spdiags(l./r,0,num_gridpts_r,num_gridpts_r); 
cos_theta — cos(theta); 
sin_theta — sin{theta); 

cos_theta_full — kron(eos_theta,ones(num_gridpts_r— 1,1)); 
sin_theta_full — kron(sin_theta,ones(num_gridpts_r— 1,1)); 

sin_theta_mat — spdiags( sin_theta ,0 , num_gridpts_theta , num_gridpts_ theta ) ; 
one_over_sin_theta — spdiags{l./sin_theta,0, ... 

num_gridpts_theta , num_gridpts_theta ) ; 



% construct divergence operator 

D— {kron(speye{num_gridpts_theta), 2*one_over_r + D_r), ... 

kron ( one_ovcr_sin_theta *D_theta*sin_theta_mat , one_ovcr_r ) } ; 



% construct gradient operator 

G — {kron(speye(num_gridpts_theta),D_r), kron{D_theta,one_over_r)}; 



% construct laplacian operators 

L — kron(speye(num_gridpts_theta),2*one_ovcr_r*D_r + D_r"2) ... 

+ kron (one_ovcr_sin_theta*D_theta*sin_theta_mat*D_th eta , one_over 



% construct surface derivative operators 

D_s — one_ovcr_sin_theta*D_theta*sin_theta_mat/r(end); 
G_s ^ D_theta/ r { end ) ; 

% construct normal derivative operator 

G_n - -kron{speye(N_t), D_r {end , : ) ) ; % d/dn = -d/dr at r = 1 



% Construct matrices to extract subsets of grid points 

% construct matrices to extract the rows corresponding to finite 
% grid points (everything except for r = infty) 

r_finite_pt_restrictor — spdiags(ones(num_gridpts_r— 1,1), 1, ... 

num_gridpts_r — 1, num_gridpts_r ) ; 
finite_pt_restrictor — kron(speye(N_t) , r_finite_pt_restrictor ); 

% construct matrix to extract the rows corresponding to interior 
% grid points {everything except for r — 1 and r — infty) 

r_interior_restrictor — spdiags(ones(num_gridpts_r— 2,1), 1, ... 

num_gridpts_r — 2, num_gridpts_r ) ; 

interior_restrictor — kron{speye(N_t) , r_interior_restrictor ); 

% construct matrix to extract the rows corresponding to r — 1 (surf 
% from a vector that already has r — infinity removed 

r_surf_restrictor — spalloc (1 , num_gridpts_r — 1 ,1); 

r_surf_restrictor(l,num_gridpts_r— 1) — 1; 

surf_restrictor — kron{speye(N_t),r_surf_restrictor); 

% construct matrix to extract the rows corresponding to r — infty 
r_inf_restrictor — spalloe(l,num_gridpts_r,l); 
r_inf_restrictor(l,l) — 1; 

inf_restrictor — kron(speye(N_t),r_inf_restrietor); 



% extract part of G operator that contributes to finite 
% using finite points 

G_f — {finite_pt_restrictor *G{ l}*finite_pt_restrictor ', 
finite_pt_restrictor *G{ 2}*finite_pt_restrictor '}; 



% split the D operators into two parts: 

% (1) contributions from finite points to finite points 
% (2) contributions from infinity to finite points 
D_f — {interior_restrictor *D{ l}*finite_pt_restrictor ', 

interior_restrictor *D{ 2 }*finite_pt_restrictor '}; 
D_inf — {interior_restrictor *D{ l}*inf_restrictor ', ... 

interior_restrictor *D{ 2} *inf_restrictor '}; 

% split the G_n operators into two parts : 
% (1) contributions from finite points 
% (2) contributions from infinity 

G_n_f — G-njt^finite.pt.restrictor '; 

G_n_inf — G_n* inf_restrictor '; 



% split the Laplacian operators into two parts 



% (1) contributions from finite points to finite points 

% (2) contributions from infinity to finite points 

L_f — intcrior_rcstrictor*L*finitc_pt_rcstrictor 

L_inf — intcrior_rcstrictor*L*inf_rcstrictor '; 

% o n t i n u at i o n loop for Newton iteration 

E = E_start; c = ones(num_gridpts,l); psi = zeros (num_gridpts,l); 
while C E <= E.final & dE > ) 

% Show progress information 

mesg — sprintf('E — %f E); disp(mcsg); 

%mm%%%%%%%%%%%%%%%%%%%%%%%%%%%mr/^^^^^ 

% compute constant terms in F = { Fl , F2 , HI , H2 ) 

wmmmmmmmmmmmmmmmmmmmmmmm^o 

Fl_const_term = c_infinity*(L_inf*ones(N_t,l)); 

F2_const_term = E*c_infinity*(— D_inf{l}*cos_theta + D_inf{2}*8in_theta) 
H2_const_term = — c _ i n f i n i t y * ( G _n _i n f * ones ( N_t , 1 ) ) ; 

% compute constant p a r 1. ol Jacobian 

DFl_Dc_const = L_f; 
DF2_Dc_const = ... 

— D_f{l}*spdiags(E*cos_thcta_full .0,num_gridpts, num_gr idpts ) ... 
+ D_f{2}*spdiags(E*sin_thcta_full ,t),num_gridpts, num_gr idpts ) ; 

DHl_Dc_const = spalloc ( N_t , num_gridpts , N_t ) ; 

DHl_Dc_COnst (: , num_gridpts_r — 1: num.gridpts.r — l:end) = . . . 

— spdiags(E*cos_theta ,0 ,N_t ,N_t); 
DH2_Dc_const ^ -G_n_f; 

% initialize loop variables using current 
% solution for c and psi 

% extract surface concentration and potential 
c_s — c ( num_gridpts_r — 1: num_gridpts_r — l:end ) ; 

phi_8 = ps i ( num_gr idp t s_r — 1: num_gr idpt s_r — 1: end ) — E*cos_theta; 

% compute zeta potential 
zeta — computeZetaPotential( ... 

V— phi_s, c_s, delta, zeta_res_tol, zeta_delta_tol, zeta_max_iters); 

% cache some common expressions 
log_c_s = log(c_s); 

s i nh_z et a_o ver _t wo = s i n h ( zet a / 2 ) ; 
cosh_zeta_over_two = cosh (zeta / 2); 

% compute surface charge density and excess neutral ion concentration 
q = — 2*sqrt(c_s).*sinh_zeta_over_two; 
w = 4*sqrt (c_s ).*( sinh (zeta / 4))."2; 

% compute initial residual 
Fl = Fl_const_term + L_f*c; 

F2 = F2_const_term + D _f { 1 } * ( c . * ( G _f { 1 } * p s i -E* c o s _t h e t a _ f u 1 1 ) ) ... 

+ D_f{2}*(c.*( G_f{2}* psi+E* sin_theta_full ) ) ; 
HI = ep s i lo n * D_s * ( q . * ( G_s * 1 o g _c _s ) + w . * ( G_S * p h i _8 ) ) ... 

— c_s . * ( G_n_f * psi + E* cos-theta ) ; 

H2 = ep8ilon*D_s*(w.*(G_s*log_c_s) + q.*(G_s*phi_s)) ... 

— G_n_f*c + H2_const_term ; 
F = [Fl; F2; HI; H2 ] ; 

res = norm(F, inf ) ; 



% Newton iteration loop 



norm_delta_soln = 1; 

count — ; 

% begin Newton iteration loop 

while (res > res_tol &&: n or m _d el t a _sol n > delta_tol count < max_iters) 

% compute Jacobian 

Dzeta_Dpsi = — l./(l+delta*sqrt(c_s).*cosh_zeta_over_two); 
Dzeta_Dc_s = — delta*sinh_zeta_over_two ./ sqrt ( c_s ) . . . 

./ (1+ delta* sqrt ( c_s ) . * cosh_zeta_over_two ) ; 
DHl_Dc_var — ( epsilon * D_s * ( ... 

+ spdiags(0.5*q./c_s.*(G_s*log_c_s) ... 

— sqrt ( c_s ).* cosh_zeta_over_two .* Dzeta_Dc_s .*(G_s*log_c_s) , . . . 
0,N_t,N_t) ... 

+ spdiags(q,0, N_t ,N_t)*G_s*spdiags(l./c_s ,0, N_t , N_t ) ... 
+ spdiags(0.5*w./c_s.*(G_s*phi_s) ... 

+ sqrt ( c_s ) .* sinh_zeta_over_two .* Dzeta_Dc_s .*(G_s*phi_s ) , . . . 
, N_t , N_t ) ) ... 

— spdiags(G_n_f*psi ,0,N_t,N_t) )*surf_restrictor; 
DHl_Dpsi_var — epsilon * D_s * ( ... 

— spdiags ( sqrt ( c_s ) . * cosh_zeta_over_two . * Dzeta_Dpsi .*(G_s*log_c_s), ... 

, N_t , N_t ) ... 
+ spdiags (w,0 , N_t , N_t )* G_s ... 

+ spdiags (sqrt ( c_s ).* sinh_zeta_over_two .* D zeta_D psi .*(G_s*phi_s), ... 
0,N_t,N_t) ) * s u r f _r e s t r i c t o r ... 

— spdiags ( c_s ,0 , N_t , N_t )* G_n_f ; 
DH2_Dc_var = epsilon * D_s * ( ... 

+ spdiags(0.5*w./c_s.*(G_8*log_c_s) ... 

+ sqrt ( c_s ) .* sinh_zeta_over_two .*Dzeta_Dc_s .*(G_s*log_c_s ) , . . . 
, N_t , N_t ) ... 

+ spdiags (w, , N_t ,N_t)*G_8*spdiags(l./c_8 ,0, N_t , N_t ) ... 
+ spdiags (0.5*q./c_s.*(G_s*phi_s) ... 

— sqrt ( c_s ) .* cosh_zeta_over_two .* D zeta _D c_S.*(G_S*phi_s), ... 
(),N_t,N_t) ) * surf_restrictor; 

DH2_Dpsi_var = epsilon * D_s * ( ... 

spdiags (sqrt ( c_s ) .* sinh_zeta_over_two .* D zeta_D psi .*(G_s*log_c_s) , . . . 
, N_t , N_t ) ... 
+ spdiags (q ,0 , N_t , N_t )* G_s ... 

— spdiags ( sqrt ( c_s ). * cosh_zeta_over_two . * Dzeta_Dpsi .*(G_s*phi_s), ... 

0,N_t,N_t) ) * surf_restrictor; 
J — [ D F 1_D c_const , spalloc(num_gridpts_interior ,num_gridpts ,0); ... 
( DF2_Dc_const . . . 

+ D_f{l}*spdiags(G_f{l}*psi ,0 ,num_gridpts ,num_gridpts) . . . 
+ D_f{2}*spdiags(G_f{2}*psi ,0,num_gridpt8,num_gridpt8) ), ... 
( D_f{l}*spdiags(c,0 ,num_gridpts , num_gridpts)*G_f{l} . . . 

+ D_f{2}*spdiags(c,0,num_gridpts,num_gridpts)*G_f{2}- ); ... 
( DHl_Dc_const + DHl_Dc_var ), D H 1 _D psi_v ar ; ... 
( DH2_Dc_const + DH2_Dc_var ), DH2_Dpsi_var ] ; 

% compute delta_soln 
dclta-soln = — J\F; 

% update solution 

c = c + delta_soln ( 1 : num_gridpts ) ; 

psi = psi + d el t a _sol n ( num_gridpt 8 +1: end ) ; 

% update residual 

'^mmmmmmmmmmmmmm tr 

% extract surface concentration and potential 



c_s = c(num_gridpts_r— l:n u m_gr i dpts_r — l:end); 

phi_s = ps i ( num_gr idpt s_r — 1: num_gridpt s_r — 1: end ) — E*cos_theta; 

% compute zeta potential 
zeta = CQmputeZctaPotential( ... 

V— phi_s , c_s , delta , zeta_re8_tol , z e t a _d e 1 1 a _t o 1 , zeta_max_iters ) ; 

% cache some common expressions 
log_c_s = log(c_s); 

8inh_zeta_over_two = sinh ( zet a / 2 ) ; 
cosh_zeta_over_two = cosh (zeta / 2); 

% compute surface charge density and excess neutral ion concentration 
q — — 2*sqrt(c_s).*sinh_zeta_over_two; 
w — 4*sqrt(c_s).*(sinh(zeta/4))."2; 

% compute residual 

Fl = F l_const_tcrm + L_f*c; 

F2 ^ F2_const_tcrm + D _f { 1 } * ( c . * ( G _f { 1 } * p s i -E* c o s _ t h e t a _ f u 11 ) ) ... 

+ D_f{2}*(c.*( G_f{2}* psi+E* sin_theta_f ull ) ) ; 
HI = eps i 1 o n * D_s * { q . * ( G_s * 1 o g _c _s ) + w . * ( G_S * p h i _s ) ) ... 

— c_s . * ( G_n_f * psi + E* cos_theta ) ; 

H2 = eps il on * D_s * (w. * ( G_8* 1 o g _c _s ) + q . * ( G_8 * phi _s ) ) ... 

— G_n_f*c + H2_const_term ; 
F = [Fl; F2; HI; H2 ] ; 

res = norm(F,inf); 

% update nor m _delta_soln , count, and residual history 
norm_delta_Soln = norm (delta_Soln ,inf); 

count — count + 1; 

% show stats 

status — [res norm_dclta_soln count] 

end % end Newton iteration loop 

% update E 
if ( E.final - E< dE) 

dE ^ E.final - E; 
end 

E = E + dE; 
end 

% Append values at infinity to results 

9mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmi> 

c — finite_pt_restrictor '*c; 
c(l:num_gridpts_r:end) — c_infinity; 
psi — finite_pt_rcstrictor '*psi; 

wmmmmmmmmmmmmmmmmmmmmmmmm^ 

% Plot results 
axis_scale = 3; 



% psi = 



potential relative to applied field 



figure(l);clf; i— i. 

min_psi = min(psi); max_psi — max (psi): ^ 

x_scale=axis_scale; y_scale=axis_scalc/2; ^ 

[rr,tt] = meshgrid(r,pi/2— theta); [ xx , yy ] = pol2cart(tt,rr); <^ 

surf (xx , yy ,reshape(psi ' ,N_r + l,N_t) '); ^sai 
hold on ; 

surf_thcta — pi/2 — [0; theta; pi]; ^ 

[surf_x,surf_y] — pol2cart(surf_theta,ones(size(surf_theta))); ^ 

surf_z = min_psi*ones(size(surf_x)); kj' 
proj.color = [180 200 220]/256; 

plots (surf_x , surf_y , surf_z , 'k ') ; fill3 (surf_x , surf_y , surf_z , proj_Color ) ; 
xlabel('x'); ylabel('z'); zlabel('\psi', 'rotation', 0); 
axis([0 x_scale — y_scale y_scale min_psi max_psi ] ) ; 

%concentration 5 

figure (2) ; clf ; *^ 
min_c = min ( c ) ; max_c — max ( c ) ; 

x_scale — axis_scalc; \'_scale — axis_scalc/2; h-j 

[rr,tt] — meshgrid{r, pi/2— theta); [ xx , yy ] — pol2cart(tt ,rr); ^^^^ 

s u r f ( XX , yy ,reshape(c',N_r + l,N_t)'); vj 

hold on; O 

surf_theta = pi/2 — [0; theta; pi]; o 
[surf_x ,surf_y] = pol2cart(surf_theta ,ones(size(surf_theta))); 
surf_z = min_c*ones(size(surf_x)); 
proj.color ^ [180 200 220]/256; 

plot3(surf_x ,surf_y ,surf_z ,'k'); fill3(surf_x ,surf_y ,surf_z ,proj_color); tS 
xlabel('x'); ylabel('z'); zlabel('c', 'rotation', 0, 'position', [ — 1.4 2 0.85]); Ot^ 
axis([0 x_scale — y_scale y_scale min_c max_c ] ) ; 

9: 
O 
O 



C.l. computeZetaPotentialO 



CD 



functionzeta = computeZetaPotential(... ^ 
Psi , c_s , delta , res.tol , de It a_z et a _t o 1 , max.iters) O 

% initialize iteration 

zeta — Psi; %uscPsiasan initial guess for zeta 

dclta_zeta = 1; norm_delta_zeta = norm (delta_zeta,inf); ^ 
res = 1; norm_res = norm (res,inf); ^ 
count = 0; 

res — zeta + 2 * d c 1 1 a * s q r t ( c _s ) . * s i n h ( z et a / 2 ) — Psi; ^■ 

CD 

% Newton iteration 

while (norm_res > res-tol & n o r m_d el t a_z et a > d e 1 1 a _z e t a _t o 1 ... 

&:count<max_iters) [3 

J = 1 + dclta*sqrt(c_s).*cosh(zeta/2); 

dclta_zcta— — res./J; I— j 
zeta = zeta + delta_zeta; O 
res = zeta + 2*delta*sqrt(c_s).*8inh(zeta/2) — Psi; 
norm_res = norm (res , inf); 

norm_delta_zet a= norm (delta_zeta ,inf); 
count = count + 1; 

end 
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